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ABSTRACT 

Self-similar solutions provide good descriptions for the gravitational collapse of 
spherical clouds or stars when the gas obeys a polytropic equation of state, p = Kp 1 
(with 7 < 4/3, and 7 = 1 corresponds to isothermal gas). We study the behaviors of 
nonradial (nonspherical) perturbations in the similarity solutions of Larson, Penston 
and Yahil, which describe the evolution of the collapsing cloud prior to core forma- 
tion. Our global stability analysis reveals the existence of unstable bar-modes (I = 2) 
when 7 < 1.09. In particular, for the collapse of isothermal spheres, which applies 
to the early stages of star formation, the I = 2 density perturbation relative to the 
background, 5p(r,t)/ p(r,t), increases as (t - t)~ - 352 (x p c (t) 0A7e , where t denotes the 
epoch of core formation, and p c (t) is the cloud central density. Thus, the isothermal 
cloud tends to evolve into an ellipsoidal shape (prolate bar or oblate disk, depending on 
initial conditions) as the collapse proceeds. This shape deformation may facilitate frag- 
mentation of the cloud. In the context of Type II supernovae, core collapse is described 
by the 7 ~ 1.3 equation of state, and our analysis indicates that there is no growing 
mode (with density perturbation) in the collapsing core before the proto-neutron star 
forms, although nonradial perturbations can grow during the subsequent accretion of 
the outer core and envelope onto the neutron star. 

We also carry out a global stability analysis for the self-similar expansion-wave solu- 
tion found by Shu, which describes the post-collapse accretion ("inside-out" collapse) of 
isothermal gas onto a protostar. We show that this solution is unstable to perturbations 
of all Z's, although the growth rates are unknown. 

Subject headings: hydrodynamics — instabilities — stars: formation - supernovae - 
ISM: clouds 



1. Introduction 

The gravitational collapse of molecular clouds leading to star formation has long been an active 
area of study. In the early stages of collapse (from p 10~ 19 g cm~ 3 to p ~ 10~ 12 g cm -3 ) the gas 
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remains approximately isothermal (at temperature ~ 10 K) due to efficient cooling by dust grains 
(see, e.g., Myhill k, Boss 1993). The gas dynamics is then specified by two dimensional parameters, 
the gravitational constant G and the isothermal sound speed a, so that the flow is expected to 
approach a self-similar form in the asymptotic limit, when the memory of initial conditions is 
"lost". Larson (1969) and Penston (1969) found a similarity solution which describes the pre- 
collapse (i.e., before the central protostar forms) evolution of the cloud, in which the gas collapses 
from rest, accelerating until it cruises at Mach number of 3.3 and the density profile reaches a 
r" 2 power law. The Larson-Penston solution contains a nonsingular homologous inner core and a 
supersonic outer envelope. A qualitatively different set of similarity solutions was found by Shu 
(1977). Of particular interest is Shu's expansion- wave solution which describes the post-collapse 
accretion of a singular isothermal gas cloud onto a protostar. In this solution, the flow starts from 
hydrostatic equilibrium (with a r" 2 density profile) and a rarefaction wave expands from the center 
and initiates the collapse (the so-called "inside-out" collapse); Inside the expansion- wave front, the 
flow eventually attains the free-fall behavior (y oc r" 1 / 2 ) at small radii, with density p oc r" 3 / 2 . 
The link between the Larson-Penston pre-collapse solution and Shu's expansion-wave solution was 
elucidated by Hunter (1977), who showed that the Larson-Penston solution can be continued to the 
post-collapse phase and that there exists an infinite (but discrete) number of pre- and post-collapse 
solutions of a different type (called "Type I"; see §2), among which the expansion- wave solution 
represents a limiting case. Figure 1 illustrates the properties of different self-similar solutions for 
the collapse and accretion of isothermal spheres. 2 

With the plethora of possible similarity solutions, it is important to know which, if any, of them 
are actually realized by collapse of isothermal clouds. One-dimensional hydrodynamical simulations, 
starting from a regular (Bonner-Ebert) sphere, generally indicate that the collapse resembles the 
Larson-Penston similarity form in the asymptotic limit (Hunter 1977; Foster & Chevalier 1993). 
This is consistent with the recent finding of Hanawa & Nakayama (1997), who showed that the 
pre-collapse Type I solutions of Hunter's (see Fig. 1) are strongly unstable against global spherical 
perturbations, and therefore are unlikely to be realized in astrophysical situations or numerical 
simulations. 

Similarity solutions have also been investigated in the context of core-collapse supernovae 
(Goldreich & Weber 1980; Yahil 1983), where the equation of state of the collapsing iron core can 
be approximated by that of a polytrope, p = Kp 1 , where K is a constant and 7 ~ 4/3. (In fact, 
the effective 7 is about 1.3 from the onset of electron capture to the neutrino trapping density, i.e., 
for 4 x 10 9 g cm" 3 <J p <J 10 12 g cm" 3 ; 7 becomes close to 4/3 when p ;> 10 12 g cm" 3 until nuclear 
density is reached.) Goldreich & Weber (1980) studied the special case of 7 = 4/3, which provides a 



2 We note that Whitworth &: Summers (1985) have found a continuum of similarity solutions by relaxing the 
analyticity condition of the flow at the transonic point; However, these generalized solutions are locally unstable 
(Hunter 1986; Ori & Piran 1988), and therefore may not be realized in astrophysical situations. We also mention 
that Boily and Lynden-Bell (1995) have constructed similarity solutions for the gravitational collapse of radiatively 
cooling gas spheres (with emissivity having a power-law dependence on density and temperature). 
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good description for the inner homologous core; They also performed a global perturbation analysis 
and showed that the inner core is stable against all radial and nonradial perturbations. Yahil (1983) 
generalized the Goldreich- Weber solution to general 7 < 4/3; this allows for a proper description 
of the outer core which collapses supersonically. Since Yahil's solution is the same as the Larson- 
Penston solution except for different values of 7, we shall often refer them as Larson- Penston- Yahil 
solutions in the remainder of this paper. 

The similarity solutions described above (in both star formation and supernova contexts) as- 
sume idealized spherical flows. A realistic gas cloud, however, contains nonradial (nonspherical) 
perturbations, 3 and it is of interest to understand the behaviors of these perturbations during 
the collapse/accretion of the cloud. In general, multi-dimensional hydro dynamical simulations are 
needed to follow the evolution of the perturbed flow, especially when the perturbations become non- 
linear. The large dynamical range involved in a collapse makes such simulations particularly chal- 
lenging (e.g., star formation ultimately involves collapse from p <, 1CT 19 g cm~ 3 to p ^ 0.1 g cm~ 3 ; 
even the initial isothermal collapse stage involves seven orders of magnitude increase in densities; 
see Truelove et al. 1997,1998 and Boss 1998 for a discussion on the numerical subtleties). An al- 
ternative, complementary approach is to carry out linear stability analysis to determine whether 
the flow is unstable to the growth of any nonradial perturbations. Since the unperturbed flow 
varies in space and time in a self-similar manner, a global analysis is needed to study perturbations 
which vary on similar scales as the unperturbed flow itself. The stability properties of the flow 
therefore depend crucially on boundary conditions at different locations of the flow. In this paper 
we perform global stability analysis for Larson-Penston- Yahil solutions (general 7) and for Shu's 
expansion-wave solution (7 = 1 only) to determine whether these similarity flows contain growing 
nonradial modes. 

While the stability properties of isothermal similarity collapse solutions (the Larson-Penston 
solution and the expansion-wave solution) are relevant to the formation of binary (and multiple) 
stars (see §5), the present study stems from our attempts to understand the origin of asymmet- 
ric supernovae and pulsar kicks (Goldreich, Lai & Sahrling 1996; see also Lai 1999). Numerical 
simulations indicate that local hydrodynamical instabilities in the collapsed stellar core (e.g., Bur- 
rows et al. 1995; Janka & Midler 1994, 1996; Herant et al. 1994), which can in principle lead to 
asymmetric matter ejection and/or asymmetric neutrino emission, are not adequate to account for 
kick velocities J> 100 km s _1 (Burrows & Hayes 1996; Janka 1998). Global asymmetric perturba- 
tions of presupernova cores may be required to produce the observed kicks. Goldreich et al. (1996) 
suggested that overstable g-modes driven by shell nuclear burning may provide seed perturbations 
which could be amplified during core collapse (see also Lai & Goldreich 2000). While the analysis 
of Goldreich & Weber (1980) shows that the inner homologous core is stable against nonradial 



3 A realistic flow may also contain a non-negligible amount of angular momentum and magnetic fields — these are 
neglected in the main text of this paper. In Appendix A we discuss the perturbative effects of rotation on Larson- 
Penston- Yahil solutions. Terebey, Shu & Cassen (1984) have considered how slow rotation affects the expansion- wave 
solution, and Galli & Shu (1993a,b) have studied the perturbative effects of magnetic fields (see also Li & Shu 1997). 
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perturbations, the situation is not so clear for the supersonically collapsing outer core where pres- 
sure plays a less important role. It is therefore important to analyse the global stability of Yahil's 
self-similar solution. Hanawa & Matsumoto (1999) have recently found a globally unstable bar 
mode in the pre-collapse Larson- Penston solution (for isothermal collapse). Our independent cal- 
culations confirm their result for 7 = 1. Since the analysis of Hanawa & Matsumoto is restricted to 
perturbations with real eigenvalues and eigenfunctions (see §3), it is not clear whether there exists 
any other growing modes, nor is it clear whether the growing bar-mode persists for general values 
of 7 (see also Hanawa k, Matsumoto 2000a) . 

The remainder of this paper is organized as follows. Section 2 summarizes the basic properties 
of the (unperturbed) Larson-Penston-Yahil similarity solution. This serves as a preparation for our 
stability analysis presented in Section 3. (For readers not interested in technical details, the main 
results are given in §3.4 and Figures 3-5.) In Section 4 we show that Shu's expansion-wave solution 
for isothermal collapse is unstable to nonradial perturbations of all angular orders. Finally, we 
discuss the astrophysical implications of our results in §5. Appendix A contains a discussion of the 
rotational and vortex modes of Larson-Penston-Yahil solutions. 



2. Spherical Larson-Penston-Yahil Self-Similar Collapse 

Here we briefly review the (pre-collapse) Larson-Penston-Yahil similarity solutions for spherical 
collapse and summarize the basic flow properties which are needed for our perturbation analysis 
(§3). 

We adopt a barotropic equation of state, where pressure p and density p are related by p = Kp 1 ', 
and K and 7 are constants. To have gravitational collapse we require 7 < 4/3. The two dimensional 
parameters of the problem are K and the Newton's constant G, from which we can construct a 
unique similarity variable 

T,= j^, R(t) = ^/2 G (l-7)/2 ( _i )( 2- 7)) (1) 

where r is the spherical radius, and the time t is measured from the epoch of core formation (i.e., the 
center formally collapses to a singularity at t = 0). Our analysis in §3 will be restricted to the pre- 
collapse solutions, so the domain of interest corresponds to t < 0. From dimensional consideration, 
we can write the dynamical (dependent) variables of the flow in self-similar forms: 



p(r,t) = Pt D( V ), Pt = G-\-t)- 2 (2) 

v(r, t) = vtV{ V ), v t = KV2 G (i-7)/2 ( _ t) i- 7 (3) 

p(r, t) = p t P(n), p t = KG-^(-t)-^ (4) 

V>(r, t) = W(r,), Vt = KG^i-tf^ (5) 

m(r,t) = m t M(ri), m t = ^/2 G (i-3 7 )/2(_ t )4-3 7 (6) 

u(r, t) = u t U( V ), u t = KG^i-tf- 2 ^. (7) 
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Here v(r, t) is the radial velocity, ip(r, t) is the gravitational potential and m(r, t) is the mass interior 
to radius r; For later purpose, we have defined the velocity stream function u such that v = Via, 
and V(rj) = dll/dr/ = U'\ here and hereafter we shall use prime (') to denote d/drj. In terms of the 
dimensionless variables, the equation of state is simply P = D 1 . The continuity equation, Euler 
equation and Poisson equation become 

WD' + DV' + 2d(i + ^) = 0, (8) 



M 



1 D^ 2 D' + WV' + ( 7 - l)y + — = 0, (9) 

M' = 4nr] 2 D, (10) 

where we have used the relation M = rj 2 ^' and have defined 

W = V + (2 - 7)7/. (11) 

Note vtW = v — (dR/dt)rj is simply the "peculiar" flow velocity with respect to the homologous 
frame. Equation (8) can be integrated out, with the help of equation (10), to give 

AnifDW = (4 - 3 7 )M. (12) 

Eliminating V' from equations (8) and (9), we obtain 

(W 2 - iDi- 1 )' 1 . (13) 



D = D 



2W 2 M 

{l-^W + 7-l)(2- 7 )r / 

7] rj z 



We see there is a sonic point at r\ = r] s , where W 2 = 7D 7 l , i.e., the (dimensionless) peculiar 
velocity W equals the sound speed (^D 1 ^ 1 ) 1 / 2 . 

Equations (10), (12) and (13) determine the spherical self-similar flow. Some properties of the 
flow are as follows. For 77 — > 0: 

D^D , V^-\n, M^^-Dov 3 ; (14) 

For rj — > 00: 

D oc r/- 2 /( 2 -T) , V oc , V/A -» = constant. (15) 

The physical solution is obtained by adjusting Dq so that the flow passes through the sonic point 
smoothly. To obtain accurate transonic solution, it is useful to analyse the behavior of the flow 
near the sonic point. The values of D s = D(r) s ), W s = W(r] s ) and M s = M(rf s ) are completely 
determined by requiring both the denominator and numerator of equation (13) to vanish at r/ s . For 
e = (rj — r] s ) /r/ s <C 1, let D = D s (l + ae) and W = W s (l + j3e). From equations (10) and (12), we 
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find M = M s [l + (2 + a + /3)e] and 2 + a + /3 
around ry s then yields 



(ry s /W s )(4 — 37). Taylor expansion of equation (13) 



(7+l)« 2 
+ 



(9 - 7 7 ) 



a + 2- (7-l)(2- 7 ) 



Vs 



(4-37)^-2 



^ + (1-7)^-4 



= 0. 



(16) 



For 7 = 1, the two roots are a = — 1 and a = rj s — 3, and the former gives the Larson- Penston 
solution. For general 7, the smaller of the two roots of (16) corresponds the Yahil-Larsen-Penston 
solution, which is the only solution with \V\ supersonic at large rj (this is called type II solution 
by Hunter 1977). The other root gives rise to a infinite (but discrete) number of solutions which 
are subsonic in \V\ at large 77 (called Type I by Hunter). We will not discuss these type I solutions 
further in this paper since they are strongly unstable against radial perturbations (Hanawa & 
Nakayama 1997). Our numerical procedure for finding the transonic solution is as follows: Guess 
Dq and r] s ; Using the boundary conditions given above, integrate equations (10) and (13) outward 
from rj = and inward from rj s to a middle point 7/ m id; Using the Newton- Rap hson scheme (Press 
et al. 1992) to vary Dq and i] s so that the two integrations match at ?? m id- 

Figure 2 gives two examples of the Larson-Penston-Yahil self-similar solutions of spherical 
collapse (for 7 = 1 and 1.3). For convenience, we list in Table 1 the key parameters of the solutions 
for different values of 7. 



3. Perturbations of Larson-Penston-Yahil Collapse Solution 

Our stability analysis relies on calculating the global linear modes of the self-similar flow. 
In general, it is not meaningful to speak of modes in flows where the unperturbed state is time 
dependent. Self-similar flows constitute an exception, since the spatial structure of the unperturbed 
flow is constant in shape, although not in scale. In this case, a mode represents a linearized 
disturbance with shape-preserving spatial structure and power-law time dependence relative to the 
unperturbed flow. The mode structure and stability depend on the feedback between boundary 
conditions at different locations of the flow. 



3.1. Perturbation Equations 

We consider flows with no net angular momentum and vorticity. 4 The fluid velocity is com- 
pletely specified by the stream function (velocity potential), i.e., v = Vu. The continuity equation, 



4 When rotation is a small perturbation, it is decoupled from the density perturbation and the potential flow. We 
discuss the rotational perturbations of self-similar flows in Appendix A. 
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Euler equation and Poisson equation for the irrotational flow can be written as 

% + V • (pV«) = 0, (17) 

^ + l (V u) 2 + h + V> = 0, (18) 

VV = 4^Gp, (19) 

where the enthalpy h = fdP/p = jKp"'^ 1 /(j — 1) for 7 7^ 1 and /i = K lnp for 7 = 1. The 
perturbed hydrodynamical equations are 

d r 1 d f 2x du\ dpd5u 2 c 



^ + ^ r J + ^ + ^ v5n = °> (20) 

—Su + + 1 Kp<~ 2 5p + # = 0, (21) 
at or 

V 2 ^ - 47rG(5p = 0, (22) 

where (5/9, <5u and (5-0 are the Eulerian perturbations of density, velocity potential and gravitational 
potential, respectively. Separating out the angular dependence in terms of spherical harmonics Y/ m , 
we write the perturbations in the form 

5 P (r,t) = (-t) s p t SD(r,)Y lin , (23) 
du(r,t) = (-t) s u t SU(ri)Y lm , (24) 
<yV(r,t) = {-t) s ip t SV(r } )Y lrn , (25) 

where pt, Vt an d u% are given by equations (2), (5) and (7). The velocity perturbation is given by 
5v(r, t) = VSu(r, t) = (-t) s v t [sV r { V ) f + 5V ± (r ] )V ± \ Y lm , (26) 

where 

^ ± =S± + -L * (27) 
at* sin if a<p 

and the dimensionless radial and tangential velocity perturbations are 

SV r ( V )=SU'( V ), SV ± ( V ) = S -^-. (28) 

In equations (23)-(26), the (unknown) power-law index s constitutes an eigenvalue of the problem. 
Since 5p(r,t)/ p(r,t) = {—t) s [5D(rj)/D{rj)]Yi m (and similarly for other variables), the value of s 
determines the global behavior of the perturbation relative to the unperturbed flow: The pertur- 
bation is globally unstable if the real part of s, Re(s), is less than zero, and it is stable if Re(s) > 0. 
Substituting (23)-(25) into equations (20)-(22), we have 

WSD' + DSU" + (2 - s + V' + 6D + (d' + ^L>) SU' - ^ — D6U = 0, (29) 
WSU' + ~fD^- 2 5D + (2 7 - 3 - s)SU + <5* = 0, (30) 
S*» = l l±V-6V - -6V + 4ttSD. (31) 

7] z 7] 
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We can eliminate SU" from (29) using (30) to obtain 



(W 2 - 7 D~<- r ) 



D 



-{3 - &y + s - 2V')6U' 



(2-s + V'+ -V) WD' 1 - 7(7 - 2)D~<- 3 D' 



SD 



+ l l±R W SU + Sty'. 



(32) 



Thus the perturbation equation is singular at the sonic point (77 = r/ s ). Equations (30)-(32) are the 
basic equations for the eigenvalue problem. 



3.2. Boundary Conditions 

To solve for the eigenvalue s, we need to know the boundary conditions. Since the unperturbed 
flow is regular at 77 — ► 0, we require the perturbations to be regular also. This gives, for 77 — ^ 0, 

SD = SD rf, SU = SU Q rf, Sty = 6V rf, (v 0) (33) 

where the three constants SDq, SUq and Styo are related by 



(t->)'- 3+2 -> 



SU - S^ . (34) 



Since the unperturbed flow is nearly static (V — > 0) at 77 — > 0, the conditions (33)-(34) are similar 
to those applied for nonradial pulsations in stars (e.g., Unno et al. 1989). 

The boundary conditions at 77 — > 00 are trickier. Let SD oc 77 s , SU oc 77^ and <5* oc T7 C for 
77 — > 00. There are four independent solutions to the fourth order systems of differential equations. 
Using the scaling relations in (15), we find that the values of a, b, c for the four solutions are 

a — 9 S 4- 2(1 — ^/) 

Solution I: o=- , b = c = a + 2 = ^; (35) 

2—7 2—7 

Solution II: a = b= - + 3 ~ 2l , c = a + 2 = * + 1 ~ 27 ; (36) 

2-7' 2-7' 2-7' y J 

Solution III: a = -I - 3 — , & = C = -(Z + 1); (37) 

2-7 

2 

Solution IV: a = /-2 , b = c = l. (38) 

2-7 

For each solution, the ratio SU/S^ and SD/S^ are uniquely determined. The general solution of 
equations (30)- (32) is a superposition of Solution I-IV. In Solution I and II, the potential perturba- 
tion S^ is produced by local density perturbation SD (thus c = a + 2); In Solution III, Sty at a large 
77 is produced by a multipole moment associated with SD at smaller 77. Solutions I-III are physically 
allowed. Solution IV, however, is not arrowed, since it corresponds to the situation where Sty at 
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a given (large) 77 is produced by density perturbation at even larger 77, and 5*$ increases without 
bound as rj — > 00. Therefore, the eigenmode at larger 77 is a linear combination of Solution I, II and 
III. Unless the Re(s) is extremely negative, i.e., for Re(s) > —(4 — 37) — (2 — 7)/, the behaviors 
of 5D and 5^ at large 77 are dominated by Solution I, while the behavior of SU is dominated by 
Solution II. Thus we have 

™ *L ^ocTf/^) (39) 

where we have used U ~ rjV oc rj^ 3 ^ 2 ^ /( 2_7 ) and ^ ~ rfD oc r/ 2 ^ 1-7 )/^ 2 " 7 ) at 77 — > 00. In practice, 
we implement the outer boundary condition at large 77 as 

5V' = -5V, c= 5 + 2(1 ~ 7) (77^00). (40) 

77 2-7 

Equation (39) indicates that when Re(s) > 0, the fractional perturbations 5D/D, SU/U and 
5^//^/ diverge as 77 increases to infinity. Thus only for globally unstable modes (Re(s) < 0) are the 
fractional perturbations finite at 77 — > 00. Whether such an unstable mode exists (for a given I and 
7) is unknown a priori. Note that equation (39) also corresponds to 

i.e., the perturbations are independent of time for 77 — > 00. 

Since the sonic point (77 = r/ s ) is a singular point of equation (32), another crucial condition 
for the perturbation analysis is that the perturbations remain regular and pass through the sonic 
point smoothly. 



3.3. Numerical Method 

Our numerical procedure for finding an eigenmode is as follows: (i) We first guess s and 
5Uo/S^>o (note that in general they are complex), and use equation (34) to find 5Dq/5^o; (In 
plotting the eigenfunctions below, we adopt the normalization such that = 1)> (h) We then 
integrate equations (30)- (32) from a small 77^ <C 1 to 773 and then from ?7 S to a large 77 ou t (we typically 
choose 77 out = 10 3 — 10 4 ); (iii) Using a Newton-Raphson scheme (Press et al. 1992), we vary the 
values of s and 5Uq/5^o until the right-hand side of equation (32) vanishes at 77,5 and condition 
(40) is satisfied at 77 ou t. Note that in step (ii), we first integrate the equations to 77 s __ = r/ s (l — e), 
where < e <C 1 (we typically choose e = 10~ 4 — 10~ 3 ), and advance the solution to 773 and to 
i] s+ = 77s (1 + e) using the derivatives evaluated at 77 s _, and then continue the integration from 
77 s+ to 77 out . We have found that this procedure works well except that for some high-order modes 
the convergence of the eigenvalue s as e decreases requires very small e. We have also tried using 
derivatives evaluated at ?7 S (and using L'Hopital's rule to calculate SD' at 773) to advance the solution 
from ?7 S _ to ?7 S _|_, but this did not lead to significant improvement. Ideally, one should not integrate 
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"into" the singular point r) s , but rather should integrate from rj s inward to a midpoint rj m i^ (< r) s ) 
and match the solution there. However, this introduces several additional unknown parameters and 
makes the multi-dimensional Newton-Raphson scheme difficult to converge in practice. 



3.4. Results 

We first note that for 1 = 1, the lowest-order mode (the one with no node in the radial 
eigenfunction) is a trivial solution; it corresponds to choosing the origin of the coordinates away 
from the center of the spherical flow. The eigenfunctions are 5D = D', 5U = U' , = The 
negative eigenvalue s = 7 — 2 should not be considered as an indication of global instability. All 
other nonradial modes are nontrivial. 



3.4.1. Unstable Modes 



For 7=1 and I = 2, we find that the lowest-order mode has a real eigenvalue, s = —0.352. 
Figure 3 depicts the eigenfunctions of the mode. Near the center (7/ — > 0), we find 5Uq/5^o = 
—0.836. The eigenfunctions are well-behaved everywhere, and go through the transonic point r] s 
smoothly. The negative eigenvalue s indicates that the bar-mode is globally unstable, with 



Sp(T, t) 

p(r, t) 
8v(r,t) 
v{r, t) 



(_f)-0.352 
(_ t) -0.352 



Y 



2m 1 



5D(rj) 
D(v) 

'sv r (n) f , g^) y ± 



v(v) 



Y, 



2m 1 



(42) 
(43) 



where p(r,t) and v(r,t) specify the unperturbed spherical flow, and 8V r {rf) = 5U'(rj), 8V±(rj) = 
5U{rj)/r). Figure 4 illustrates the growth of the density perturbation as the collapse proceeds. The 
fractional perturbation grows as (— t) -0 - 352 oc p c (i)°- 176 , where p c (t) = p(0,t) is the central density 
of the cloud. The growing bar-mode corresponds to the deformation of the collapsing cloud toward 
an ellipsoidal shape. Depending on the initial perturbations, the deformed cloud may take the form 
of an oblate disk or a prolate bar (see also Hanawa & Matsumoto 1999). 

As 7 increases, the mode tends to be stablized by the effect of pressure. Figure 5 depicts the 
variation of s = sq for the lowest-order bar-mode (I = 2) as a function of 7. We find that s increases 
with increasing 7, and the mode is unstable (with negative s) only for 7 < 1.09. 5 Figure 6 gives a 
few examples of the mode eigenfunctions for several different values of 7. 



5 Similar result is also obtained by Hanawa & Matsumoto (2000a) in a different analysis. The author thanks the 
referee, T. Hanawa, for pointing out this paper. 



3.4-2. Stable Modes 



We have searched numerically for other unstable modes [with Re(s) < 0] for 1 < 7 < 4/3 and 
I = 1, 2, 3, • • •. Our search covers the domain —5 <^ Im(s) <^ 5. However, except for those discussed 
in §3.4.1, all modes we have found are stable [with Re(s) > 0]. As an example, the dashed curve 
in Fig. 6 shows the eigenfunction of a high-order I = 2 mode (for 7 = 1), with s± = 0.23 + 0.261 
Note that as 7 increases, the ordering of the modes can change. This is seen from Figure 5: For 
7 ^ 1.11 we have sq <Re(si), but for 7 ;> 1.11 we find sq >Re(si). We have not explored the 
spectrum of high-order modes in detail, since these modes are all stable. Moreover, as the fractional 
perturbations associated with the stable modes diverge in the r? — > 00 limit (see eq. [39]), these 
modes are only formally well-defined, but are of no physical importance. 



4. Perturbations of "Inside-Out" Collapse of Isothermal Cloud 

In this section we present our perturbation analysis of Shu's expansion-wave solution which 
describes the "inside-out" collapse of a isothermal gas cloud. The equation of state is p = Kp = pa 2 , 
where a is the sound speed. 



4.1. Spherical Inside-Out Collapse: Shu's Expansion- Wave Solution 



The expansion-wave solution describes the post-collapse (t > 0) evolution of the flow. The 
similarity variable is defined as 

r] = —. (44) 



at 



The flow variables can be written in self-similar forms as in equations (2)-(7), except that in 
pt, vt, pt, • •• we have to replace (— t) by t and set 7 = 1, i.e., 



1 



aH 



Pt AirGt 2 ' 



v t = a, p t = ^ , ipt = a , m t = — , 



ut = a t. 



(45) 



(Note that to follow Shu's convention, we have included the factor 4tt in pt.) In terms of the 
similarity variables, the continuity equation, Euler equation, and Poisson equation are 



(V - rj)D' + DV' + 2d(^--1^= 0, 

D 7] z 

M' = r] 2 D. 

These equations can be rearranged into the standard form as given by Shu: 



[(v-v) 2 -i}^ = (v-v) 



D- 2(1-- 



(46) 

(47) 
(48) 



(49) 
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[(V- V ) 2 -l] V'=(tj-V) 



D(v -V)-- 
77 J 



(50) 



and M = rf(r)- V)D. 



Some properties of the expansion-wave solution are as follows. For 77 > 1, the solution describes 
a static isothermal sphere, with V(i]) = and D(r}) = 2/rj 2 . The surface 77 = 1 is the rarefaction 
(expansion) wave front. For r/ — > 0, the solution describes a free-fall, with M — > Mq = 0.975, 
V - (2M /r]) 1/2 , and L> -» (M /2r/ 3 ) 1/2 . While L> and F are continuous at 77 = 1, £>' and V' 
are not: 

V'(l+)=0, £>'(l+) = -4; V'(l-) = 1, £>'(l-) = -2. (51) 
(The notation 77 = 1+ means that 77 — ^ 1 from above, and 17 = 1— means 77 — > 1 from below.) 



4.2. Perturbation Equations 

As in equations (23)- (25), we consider perturbations of the form 

Sp(r,t)=t s p t SD( V )Y lm , 
Su(r,t)=t s u t SU( V )Y lrn , 
S^(r,t)=t s ^ t S^( V )Y lm . 



(52) 
(53) 
(54) 



Since 5p(r,t)/p(r,t) = t s [SD(r])/D(rj)] Y\ m (and similarly for other variables), the power-law index 
s specifies the evolution of the perturbation relative to the background: The flow is unstable if 
Re(s) > and stable if Re(s) < 0. Substituting (52)-(54) into the perturbation equations (20)- 
(22), we obtain 

(V - j])8D' + DSU" +(-2 + s + V' + SD + [d' + ^) Su ' ~ ^ ^ D6U = °> ( 55 ) 

(56) 



SD 



(V - n)SU' + — + (1 + s)5U + 6V = 0, 

Sy» + - ^±^<5tf -SD = 0. 
77 77^ 

We can use equation (56) to eliminate SU" in equation (55) and obtain 



(57) 



(V-rif-1 ^ - (2V' + s)5U' + 



D 

1(1 + 1) 

7] 2 



(-2 + S + V+ lv}{y-ri) + 



iy 
15 



5D_ 

IT 



(V - r])5U - <5*' = 0. 



(58) 



Thus, the expansion-wave front (77 = 1) is a singular point of the perturbation equation. Also, 
equation (57) can be written in the integral form: 



= -v l P(v) - ^r- 



(59) 
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where 

P(n\ = , , , . 

21 + 1 J v ' w ' " 2/ + 1 



1 Z -00 1 
P(n) = / r,' l - l 5D{r,')dr}', P' = rj^SD, (60) 



Q(V) = ^J\ ,l+2 SD(r,')d V ', Q' = -L-r, l+2 5D. (61) 



4.3. Series Solution for rj > 1 

For rj > 1, we have F = and D = 2/rj 2 , the perturbation equations can be solved in Frobenius 
series. We consider the solution which satisfies 5D/D — ► 0, SU — > 0, and 5^> oc ?7 _ ' _1 — > for 77 — > 00 
(i.e., is given by the decreasing solution of the Laplace equation); 6 The last condition implies 
that Q approaches a constant as rj — > 00. Thus we can write 

00 

Q(v) = ^2Q2nV- 2n - (62) 

n=0 

Equation (61) then gives 

00 

5D(rj) = d2nV~ 2n ~ l ~ 3 , d 2n = -2n(2Z + 1) q 2n . (63) 

n=0 

Using equations (59), (60) and requiring P — > as rj — > 00, we have 

00 2/ + 1 

= 5>n?T 2n -'-\ V>2n = ~ 2n+ ^ + 1 g2n. (64) 

n=0 

Substituting (63) and (64) into (56) yields 

, TT( x V 2 - -2n-l-i (2l + l){2n 2 +2ln + n+l) 

6U(rj) = \ u 2n = ( L + 2 ;Vl)(2n + / + 2 + S V 2 - (65) 

n=0 v 7 v ' 

Finally, using equation (58), we obtain the recurrence relation: 

(2n + 3 + / + s)d 2 „+ 2 = (2n + I + l)d 2n + 2 [s(2n + / + 1) + /(/ + 1)] u 2 „ 

+2(2n + / + l)V>2n (n = 0,l,2,---) (66) 

With this recurrence relation, the complete solution for rj > 1 can be obtained. Note that for 
rj — ► oo, the asymptotic scalings of the perturbations are 

SV^-^, SU^- — — — - — ttt , 5D (x-}—. (67) 

rj l+1 (l + 2 + s)rj l+i rj l + 5 V ; 



6 The fourth order system of differential equations allows for four independent solutions, but this solution (which 
must exist for any physical situation) alone is adequate for our stability analysis (§4.4). 
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4.4. Instability 

Here we use the series solution of §4.3 and the boundary condition at the expansion-wave front 
(77 = 1) to show that Shu's solution is unstable. As equation (58) indicates, the expansion-wave 
front is a singular point of the perturbation equation. A natural (necessary) boundary condition 
at rj = 1 is that the perturbation is finite (although SD and SV r = 5U' can be discontinuous across 
77 = 1; see below). 

We can examine the behavior of the perturbation at rj — > 1+ using the series solution of §4.3. 
From the recurrence relation (66) we find, for n — > 00, 

d 2n+ 2 , , l + S (6g) 



d 2n n 
u 2n +2 2 + 3 

u 2n n 

V>2n+2 _^ 1 _ 3 + g 

1p2n U 



(69) 
(70) 



Thus in order for 5D to be finite at rj — > 1+, we require Re(s) > (e.g., Mathews & Walker 
1970). One can similarly show that in order for 5V r = SU' to be finite at 77 — > 1+, we require 
Re(s) > 0. Thus, any perturbations which are well-behaved at the expansion-wave front must be 
globally unstable. 

A possible caveat in the analysis given above is that in the presence of flow perturbations, 
the rarefaction front is also perturbed, and 5D(r] — > 1+) does not give the density perturbation 
at the perturbed expansion-wave front; one might therefore be concerned that the divergence of 
5D(rj — ► 1+) is a result of an improper definition of 5D. To address this problem, we define a 
stretched radial coordinate via 

£(77) = V (1 + AY lm t s ) , (71) 

where A is a constant (to be determined). The perturbed rarefaction front is located at £(7/ = 1). 
Since D(rj) + 5D( V )Y lm t s = D[£(rf)] + 5D[^ V )}Y lm t s , we have 

5D[£(v)] =5D(r))-r)D\n)A. (72) 

Similarly, SU'[^(r))] = 5U'{rj) — rjV'(rj)A. Since 5D[£(r))] and 8U'[£(r))] must be continuous across 
the rarefaction front, and since D' and V are discontinuous at rj = 1 (see eq. [51]), we infer that 
5D(rj) and 5V r (rj) are discontinuous at rj = 1. Evaluating equation (58) at 77 = 1+ and rj = 1—, we 
find 

5U'{1+). (73) 



8+1 

Using equation (72) we obtain 



5D[i(r, = 1+)] = 6D(r) = 1+) - —^—SU\f, = 1+). (74) 

(s + 1) 
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Using the series solution of §4.3, we can easily show that 8D[£(r] = 1+)] diverges unless Re(s) > 0. 

Another concern one might have is that the divergence of 5D and 5V r at rj = 1+ discussed above 
simply indicates that the series expansion breaks down at rj = 1 rather than the actual divergence 
of the function 5D and SV r . To address this issue, we show in Figure 7 several examples of the 
absolute value of the density perturbation 5D at small {rj — 1) for several different values of s. The 
function 5D is calculated using the series expansion given in §4.3 (normalized by setting qo = 1). 
We see that, in accordance with our discussion above, when Re(s) < 0, the density perturbation 
5D diverges as 77 — > 1+. Indeed, an analysis of the perturbation equations near 77 = 1+ shows that 
for < x = 77 — 1 < 1 the perturbations have the following behavior: 

SD = C x s [1 + 0(x)\ + Ci [1 + 0(x)\ , (75) 

5U = 2(7+T) xS+1 [1 + + ° 2 [1 + ' (76) 
= (7TT^+2) xS+2 [1 + ° (a;)] + ° 3 [1 + ' (77) 

where Co, Ci, C2, C3 are constants. This clearly shows that 5D(r] = 1+) diverges for Re(s) < 
- We could have deduced this result simply by examing the perturbation equations near 77 = 1+, 
except that without the series solution discussed in §4.3 we would not know whether Co = is a 
possibility. The numerical results (based on the series solution) depicted in Figure 7 agree with 
(75)-(77) and Co / 0, i.e., the boundary condition at 77 — ► 00 requires Co 7^ 0. It is this global 
consideration of the perturbations at 77 — ► 00 and at 77 — > 1+ that forces us to conclude that the 
expansion-wave solution is unstable to perturbations of all Vs. 

Note that our analysis above indicates Re(s) > 0, but we have not solved for s. (The actual 
values of s depend on the flow at 77 < 1 and the boundary conditions at 77 — > 0.) Thus the growth 
rates of the instabilities are unknown at present. 

5. Discussion 

Early studies by Hunter (1962) and by Lin, Mestel & Shu (1965) demonstrated that uniform, 
pressure-free gas clouds undergoing gravitational collapse are unstable toward fragmentation and 
shape deformation, with perturbations growing asymptotically as 5p(r, t)/p{t) oc (to—t)^ 1 oc pit) 1 / 2 
in the linear regime, where to denotes the epoch of complete collapse, and p(t) is the unperturbed 
uniform density. However, the presence of even a small initial central concentration and pressure 
forces significantly alters the evolution of the cloud. If the gas pressure is simply related to the 
density by a power-law, p = Kp 1 (polytropic equation of state), the flow asymptotically approaches 
the similarity solutions found by Larson (1969), Penston (1969) (for isothermal gas 7 = 1), by 
Goldreich Sz Weber (1980) (for 7 = 4/3), and by Yahil (1983) (for general 7). Since the local Jeans 
length is of the same order as the length scale at which the flow varies, a global analysis is needed 
to determine the stability properties of the collapsing cloud. The result (§3) presented in this paper 
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(see also Hanawa & Matsumoto 1999) shows that for sufficiently soft equation of state (7 < 1.09), 
the Larson-Penston-Yahil similarity flow is unstable against bar-mode perturbations, such that 
5p(r,t)/p(r,t) oc (to — t) s Y2 m (6, 4>) with s < (s = —0.352 for 7 = 1 and s increases to zero as 
7 increases to 1.09, see Fig. 5), where to denotes the epoch of core formation. Since the central 
density increases as p c (t) oc (i — t)~ 2 , the growth of perturbation, Sp(r,t)/p(r,t) oc p c (t)~ s / 2 , is 
slow (e.g., for isothermal collapse, 5p/p increases by a factor of 1.5 when p c increases by a factor 
of 10). Such a slow growth (compared with the Sp/p oc p 1 / 2 behavior for the collapse of uniform, 
pressure-less gas) is a result of the stablizing influence of pressure, despite the large Mach number 
(about 3) achieved in the outer region of the cloud. 

Our stability analysis applies to the pre-collapse stage (prior to core formation) of the Larson- 
Penston-Yahil solutions. After the central core forms, the outer core and envelope accrete onto it 
(see Fig. 1). The gas approaches free-fall as r — ► 0, and the Mach number becomes much greater 
than unity. In this (accretion) stage, nonradial perturbations (of all scales) grow kinematically as 
Sp/p oc r -1 / 2 oc p 1 / 3 , where r(t) is the radius of a fluid element and p(t) oc r~ 3 / 2 its comoving 
density (Lai & Goldreich 2000). Although the fluid element is free- falling, the perturbation grows 
more slowly compared with the case of uniform pressure-less collapse because the steep velocity 
gradients provide a stablizing influence on the flow. 

The global bar-mode instability for isothermal collapse may have important implications for 
star formation, particularly in connection with the formation of binary (and multiple) stars (see 
also Hanawa & Matsumoto 1999; Matsumoto & Hanawa 1999). Fragmentation is unlikely to occur 
in a globally spherical collapse because small condensations do not contract fast enough to separate 
out from the converging bulk flow. Angular momentum (or magnetic field) can obviously make the 
cloud nonspherical, and thus facilitate fragmentation (e.g., Burkert & Bodenheimer 1996; Burkert, 
Bate & Bodenheimer 1997; Truelove et al. 1997,1998; Boss 1998). Observations suggest that 
many of the molecular cloud cores (with mass of order a few M Q and size 0.1 pc) have elongated 
shapes (Myers et al. 1991) and slow rotation rates (with the ratio of rotational to gravitational 
energies of order 0.02; Goodman et al. 1993), implying that rotation is probably not a crucial 
factor in driving fragmentation on scales greater than 200 AU. Our result on the growth of bar- 
mode perturbations {5p oc Yim) indicates that, even without net angular momentum, the collapsing 
cloud tends to deform into an ellipsoidal shape (oblate disk or prolate bar, depending on which m- 
mode perturbation is dominant initially). Fragmentation is more likely to occur for such deformed 
configurations (e.g., Bonnel 1999; Matsumoto &; Hanawa 1999). 

In the context of core-collapse supernovae, our result shows that the homologous inner core and 
the supersonic outer core are globally stable against nonradial perturbations prior to core bounce 
at nuclear density and the formation of the proto-neutron star. However, during the subsequent 
accretion of the outer core (involving 15% of the core mass) and envelope onto the proto-neutron 
star, nonspherical perturbations can grow according to 5p/p oc r -1 / 2 or even bpj p oc r _1 (Lai & 
Goldreich 2000). The asymmetric density perturbations seeded in the presupernova star, especially 
those in the outer region of the iron core, are therefore amplified during collapse. The enhanced 
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asymmetric density perturbation may lead to asymmetric shock propagation and breakout, which 
then give rise to asymmetry in the explosion and a kick velocity to the neutron star (Goldreich et 
al. 1996; Burrows & Hayes 1996). 

Our stability analysis (§4) shows that Shu's expansion-wave solution is globally unstable to 
perturbations of all Us, although the growth rates are unknown at present. The implication of this 
result is not entirely clear. It is well-known that a static singular isothermal sphere is highly unstable 
to radial perturbations (A truncated Bonner-Ebert isothermal sphere is unstable when the range of 
density from the center to the surface is greater than 14.04; see Bonner 1956, Hunter 1977). Earlier 
one-dimensional numerical simulations have already shown that a collapsing isothermal cloud does 
not approach the expansion- wave solution (Hunter 1977; Foster Sz Chevalier 1993). Our stability 
analysis corroborates this result, and indicates that the expansion-wave solution cannot be realised 
in a pure hydrodynamical situation. 

Magnetic fields play an important role in the current paradigm for forming low- mass stars (e.g., 
Shu, Adams & Lizano 1987; Shu et al. 1999; Mouschovias & Ciolek 1999). Ambipolar diffusion of 
magnetic fields drives the quasi-static contraction of the molecular cloud core with growing central 
concentration such that the core asymptotically approaches the state of a singular isothermal sphere. 
When the flux-to-mass ratio drops below certain critical value, a runaway "inside-out" collapse 
ensues, and it is thought that this collapse is well described by the expansion-wave solution (Shu et 
al. 1999). In reality, there is probably no sharp distinction between the quasi-static contraction and 
dynamical collapse (e.g., Safier, McKee & Stahler 1997; Li 1998), and a real singular isothermal 
sphere can never be reached. Our global stability analysis of the expansion-wave solution (§4) does 
not depend on the mathematical singularity of the solution at r = 0, but depends on the existence 
of a well-defined rarefaction front and a static isothermal density profile outside the front in the 
solution. It in not clear whether our idealized hydrodynamical stability analysis can be applied to 
more realistic situations with (even sub-dominant) magnetic fields (see Galli & Shu 1993a, b and Li 
& Shu 1997 for the effects of magnetic field on self-similar "inside-out" collapse). 

This work was started in 1995 when I was a postdoc in theoretical astrophysics at Caltech 
(support from a Richard C. Tolman fellowship is gratefully acknowledged). I thank Peter Goldreich 
for initially suggesting this problem in the context of core-collapse supernovae and for many valuable 
discussions. I also thank Frank Shu and the referee, T. Hanawa, for useful comments on this paper. 
This work is supported in part by NASA grants NAG 5-8356 and NAG 5-8484, and by a research 
fellowship from the Alfred P. Sloan foundation. 

A. Rotational Perturbations in Larson-Penston-Yahil Solutions 

A general velocity perturbation can be written as 

Sv(r,t) =Sv r (r,t)Y lm f + 5v ± (r,t)V L Y lm + V ± x [5v TOt (r,t)Y lm f] . (Al) 
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Using Euler equation, we obtain 

j t (r6v Iot ) =0, (A2) 

Jt 6VT = -lfr 5VT ' (A3) 

where Svt = 5v r — d(r5v±)/dr (see Lai k Goldreich 2000). The potential flow discussed in the 
main text corresponds to 5u = r5v_\_, 5vt = and 5v rot = 0. Note that 8v rot is decoupled from the 
potential flow. 

Writing 6v TO t in the self-similar form, 5v rot (r,t) = vt(—t) s SV ro t(rj), equation (A2) becomes 

W( V SV Iot )' + (s - 3 + 2 7 )(r / <5V rot ) = 0, (A4) 

where W = V + (2 — 7)77. Since F oc rj as 77 — > 0, it is most natural to require SV rot oc r\ at r\ — > 
(corresponding to a uniform "rotation"). Equation (A4) then gives s = —1/3, independent of 7. 
This is a growing mode which describes the spin-up of a rotating cloud during gravitational collapse 
(see Hanawa & Nakayama 1997; Matsumoto & Hanawa 1999). The "angular frequency" increases as 
$v TO t/r oc (— t) -4 / 3 <5T4ot/77) and the velocity perturbation increases as 5v Iot /v cc {-ty 1 / 3 oc pc(t) 1 / 6 . 

Similarly, writing Svt as 5vT(r,t) = vt(—t) s 5Vr(ri), equation (A3) becomes 

WSVt + (V + 7 - 1 - s)5F T = 0. (A5) 

For 77 — > 0, we have <5V^. oc r/" 1 , <5Vj_ oc r/' -1 , but 5Vr oc r/' +1 . Equation (A5) then gives s = 
(4/3 — 7)/ — 1/3. This is the growing "votex" mode discussed by Hanawa & Matsumoto (2000b). 
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Table 1. Parameters for Pre-collapse Larson-Penston-Yahil Solutions 





7 




Vs 




1 


.00 


0.13256 


2.34113 


3.271 


1 


.05 


0.18299 


2.33723 


2.802 


1 


.10 


0.25908 


2.33252 


2.506 


1 


.15 


0.37769 


2.32991 


2.338 


1 


.20 


0.57228 


2.33277 


2.290 


1 


.25 


0.92455 


2.34606 


2.407 


1 


.30 


1.75375 


2.38512 


2.944 


1 


.31 


2.10358 


2.40203 


3.216 


1 


.32 


2.66174 


2.42909 


3.695 


1 


.33 


3.99500 


2.49457 


5.113 



Note. — 7 is the polytropic in- 
dex, Dq = D(j] = 0), Tj s is the 
sonic point, where W = A = 
( 7 D7-i)i/2, and M 0O = (\V\/A) 0O 
is the Mach number at rj — > oo. 
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Fig. 1. — A schematic diagram showing the properties of different self-similar solutions describing 
the collapse and accretion of isothermal gas clouds. The similarity variable is r/ = r/(—at) for 
pre-collapse solutions (t < 0) and rj = r/{at) for post-collapse (accretion) solutions (t > 0), where 
a is the sound speed, and t = corresponds to the epoch when the central core collapses to form 
a protostar. The vertical axis gives the dimensionless inflow flow velocity —V = —v/a. Shu's 
expansion- wave solution (the solid curve that terminates at rj = 1, the rarefaction wave front) 
describes post-collapse accretion, while the other solutions have a pre-collapse phase and a post- 
collapse phase which are connected at rj — ► oo (or t = 0). All post-collapse flows approach free-fall 
V oc r/ -1 / 2 as rj — > 0. At rj — > oo, the Larson-Penston solution has Mach number of 3.3. The dashed 
curves give an examples of the infinite (but discrete) number of type I solutions found by Hunter 
(1977). (Note that the pre-collapse type I solutions contain regions with both positive and negative 
v.) The expansion-wave solution is the limiting case (V — > at i] — > oo) of the post-collapse Type 
I solutions. Note that all pre-collapse solutions have V — > —2i]/3 as rj — > 0. 
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Fig. 2. — The Larson-Penston-Yahil similarity solutions are shown for 7 = 1 and 7 = 1.3. The solid 
curves give the dimensionless flow velocity (— V), and the dashed curves give the density profile D. 
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Fig. 3. — Eigenfunctions of the lowest-order bar-mode (I = 2) for 7 = 1 (isothermal collapse), with 
the eigenvalue s = —0.352. The upper panel shows the fractional density perturbation SD/D, the 
middle panel shows the radial and tangential velocity perturbations, and the lower panel shows 
the potential perturbation. The dotted vertical line denotes the transonic point rf s = 2.341. The 
similarity variable is 77 = r/(—at), where a is the (isothermal) sound speed. 
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Fig. 4. — Evolution of the I = 2 density perturbation during the collapse of an isothermal 
cloud (7 = 1). The angular dependence, i^m, has be suppressed. Note that 5p(r,t)/ p(r,t) oc 
(—t)~°- 352 5D(ri)/D(r]), with 77 = r/(—at); T is a fiducial time, and a is the sound speed. The differ- 
ent curves correspond to different times. The center of the cloud reaches singularity as t approaches 
zero. 
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Fig. 5. — The eigenvalue s of the bar-mode (I = 2) as a function of the polytropic index 7. The 
filled circles correspond to the lowest-order mode (with no radial node in the eigenfunction) , with 
s = sq real; The open circles correspond to Re(si) of a higher-order mode, with s = s\ complex, 
the dashed curve gives Im(si) of the same mode. Note that the lowest-order bar mode is globally- 
unstable for 7 < 1.09. The s = si mode has one radial node for 7 < 1.15, and it crosses the 
zero-node mode (s = so) at 7 ~ 1.11. For 7 > 1.11, the s = si mode is the mode with the lowest 
Re(si). 
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Fig. 6. — The eigenfunctions 5D/D of several bar-modes (I = 2) for different 7. The solid curves 
correspond to the lowest-order bar-mode for 7 = 1 (with s = —0.352, unstable), 7 = 1.08 (with 
s = —0.038, unstable), and 7 = 1.1 (with s = 0.106, stable). The dashed curve gives Re(5D/D) 
for the s = s\ = 0.23 + 0.26i mode (see Fig. 5) with 7 = 1. 
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Fig. 7. — Behavior of the absolute value of the density perturbation 5D just outside the expansion- 
wave front (77 = 1) in the expansion-wave solution. The solid curves are for I = 2 and the dashed 
curve for I = 1. The values of s are labeled for each curve. Note that when Re(s) < 0, the 
perturbation \5D\ — > 00 as 77 — > 1. 



